#ifndef __vertex_fe_hh__
#define __vertex_fe_hh__

#include "nodal_fe.hh"
#include "fem_traits.hh"

/// Class for constant FE on interval
template<typename T=double>
class VertexFiniteElement : public NodalFiniteElement<0, fastfem::Vertex, 0, fastfem::Lagrange, T> 
{
public:
   VertexFiniteElement() : NodalFiniteElement<0, fastfem::Vertex, 0, fastfem::Lagrange, T> () {}
   virtual ~VertexFiniteElement() {}

  /** virtual function which evaluates the values of all
       shape functions at a given point ip */
  virtual T evaluate_basis (  unsigned int , const Point<SPATIAL_DIMENSION, T>& )
  {
    return 1.0;
  }
  
  /** virtual function which evaluates the values of all
       shape functions at all point ips */
  virtual void evaluate_basis_all ( const PointList<SPATIAL_DIMENSION, T>& ips, std::vector<T>& basis_vals )
  {
    basis_vals.resize(ips.n_points, 1.0);
  }
  
  virtual void evaluate_gradient (  unsigned int , const Point<SPATIAL_DIMENSION, T>&, Point<SPATIAL_DIMENSION, T>& grad)
  {
    grad.reset();  // all gradients are zero
  }
  
  virtual void evaluate_gradient_all ( const PointList<SPATIAL_DIMENSION, T>& ips, std::vector<Point<SPATIAL_DIMENSION, T> >& grad_vals )
  {
    grad_vals.resize(ips.n_points, Point<SPATIAL_DIMENSION, T>()); // all gradients are zero
  }
  
  virtual void evaluate_hessian ( unsigned int , const Point<SPATIAL_DIMENSION, T>&, PointList<SPATIAL_DIMENSION, T>& hess)
  {
    hess.reset() ;
  }

  virtual void evaluate_hessian_all ( const PointList<SPATIAL_DIMENSION, T>& ips, std::vector<PointList<SPATIAL_DIMENSION, T> >& hess_vals )
  {
    hess_vals.resize(ips.n_points, PointList<SPATIAL_DIMENSION, T>()); // all gradients are zero
  }

  void reference_element_coords(std::vector<Point<SPATIAL_DIMENSION, T> >& nodes)
  {
    nodes.resize(1);
    nodes[0][0] = 0.0;  // origin
  }
  
  virtual int test()
  {
    std::vector<Point<SPATIAL_DIMENSION, T> > nodes ;
    this->reference_element_coords(nodes);
    
    cout << "\n -- VertexFiniteElement --";
    cout << "\n   Number of reference element nodes:= " << nodes.size() ;
    cout << "\n        Reference element coordinates = [" << nodes[0][0]  << ", " << nodes[0][1]  << ", " << nodes[0][2] << "]" << endl ;
    
    return 0 ;
  }
  
};


template<typename T>
class FEMTraits<VertexFiniteElement<T> >
{
  
  enum { FE_Dimension = 0 };
  enum { FE_Family = fastfem::Lagrange };
  enum { N_Nodes = 1 };

  typedef T BaseType;
  typedef VertexFiniteElement<T> FEType;
  typedef NodalFiniteElement<0, fastfem::Vertex, 0, fastfem::Lagrange, T>  ParentFEType;

};


#endif // #ifndef __vertex_fe_hh__
